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Fluctuating hydrodynamics has been successfully combined with several computational 
methods to rapidly compute the correlated random velocities of Brownian particles. In the 
overdamped limit where both particle and fluid inertia are ignored, one must also account 
for a Brownian drift term in order to successfully update the particle positions. In this pa¬ 
per, we present an efficient computational method for the dynamic simulation of Brownian 
suspensions with fluctuating hydrodynamics that handles both computations and provides 
a similar approximation as Stokesian Dynamics for dilute and semidilute suspensions. This 
advancement relies on combining the fluctuating force-coupling method (FCM) with a new 
midpoint time-integration scheme we refer to as the drifter-corrector (DC). The DC resolves 
the drift term for fluctuating hydrodynamics-based methods at a minimal computational 
cost when constraints are imposed on the fluid flow to obtain the stresslet corrections to the 
particle hydrodynamic interactions. With the DC, this constraint need only be imposed 
once per time step, reducing the simulation cost to nearly that of a completely deterministic 
simulation. By performing a series of simulations, we show that the DC with fluctuating 
FCM is an effective and versatile approach as it reproduces both the equilibrium distribu¬ 
tion and the evolution of particulate suspensions in periodic as well as bounded domains. 
In addition, we demonstrate that fluctuating FCM coupled with the DC provides an effi¬ 
cient and accurate method for large-scale dynamic simulation of colloidal dispersions and 
the study of processes such as colloidal gelation. 

Keywords: Brownian motion, fluctuating hydrodynamics, Stokesian dynamics. Brownian 
Dynamics, colloidal suspensions, simulation 


I. INTRODUCTION 

Brownian motion is the random motion ex¬ 
hibited by micron and sub-micron particles im¬ 
mersed in liquid. It is a fundamental mechanism 
for material and chemical transport in micron- 
scale physical and biological system^, and can 
play a key role in determining the mechani¬ 
cal res ponse of colloidal suspensions to applied 
stressej^l ^bi | 24 |_ Brownian motion is also known 
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to affect the aggre gation and self-assembly of 
interacting particleP®^^, a fundamental pro¬ 
cess important in many engineering applications 
that utilize colloidal p articl es to tune rheolog¬ 
ical properties of flui dj^^*^ and construct new 
materials and deviceP^^^^^^. 

While Brownian motion is clearly important 
in these situations, the careful study and quan¬ 
tification of its role using simulation remains 
a computational challenge due to the intimate 
link between the random motion of the par¬ 
ticles and their many-body hydrodynamic in¬ 
teractions. For suspensions of moderate con¬ 
centration, this requires resolving the hydrody¬ 
namic interactions beyond those provided by 
point forces. It involves obtaining the second- 
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order corrections to the hydrodynamic interac¬ 
tions that are provided by the stresslets - the 
symmetric force moments on the particles - as 
well as the rotlets due to torques applied to the 
particles. The stresslets are necessary to quan¬ 
tify suspension rheology and accurately resolve 
particle dynamics in linear flow fields. 

The inclusion of stresslets in Brownian sim¬ 
ulations presents several computational chal¬ 
lenges and has been addressed directly in Stoke¬ 
sian Dynamics (SD), and in accelerated SD 
through approximation. The primary contribu¬ 
tion of this work is a complete method to per¬ 
form dynamic Brownian simulations of dilute 
and semidilute suspensions that includes the ef¬ 
fects of stresslets exactly using fluctuating hy¬ 
drodynamics. We do this in a way that incurs 
only a minimal computational cost over a de¬ 
terministic simulation of the same system, and 
in a framework that we show can be extended 
beyond periodic boundary conditions. 

In order to achieve the correct equilib¬ 
rium distribution, the statistics of the random 
particle motion must satisfy the fluctuation- 
dissipation theorenP^ which states that the ran¬ 
dom particle velocity correlations must be pro¬ 
portional to the hydrodynamic mobility matrix. 
In simulation, this would require the square 
root of the mobility matrix to be found at each 
time step to compute the correct particle ve¬ 
locities. As a result, the inclusion of Brownian 
motion has often limited simulations to having 
very small particle numbers, or ignoring com¬ 
pletely hydrodynamic interactions between the 
particles. Though several methods have been 
introduced t o acce lerate the matrix square root 
computatiorP^^^, recent studies have shown 
this computation can be avoided altogether us¬ 
ing fluctuating hydrodynamics. 

Fluctuating hydrodynamics involves generat¬ 
ing random fluid flows by including a white- 
noise fluctuating stress in the equations of 
fluid motion. Introduced in the first edi¬ 
tion of Landau and LifshitJ^, its effective¬ 
ness for yielding the correct random mo¬ 
tion of particles was demonstrated in a num¬ 
ber of theoretical studies in the 1960s and 
1970 j 9 | iu | i 4 | 25 | 32 jllIMI 7 T]^ As a result of this 


fundamental work, fluctuating hydrodynam¬ 
ics has found success in numerical simulations 
of micron-scale fluid-structure interactions in 
methods such as Lattice-BoltzmaniPHUl^ hy¬ 
brid Eu lerian -Lagrangian approaches for point 
particle^SSmi^ distributed Lagrange-multiplier 
methocP^, finite-element simulatioiP^, the 
stochastic and f luctuati ng immersed-boundary 
methods and the fluctuating 

force-coupling method (FCMj^. 

IBM and FCM are similar in that they both 
use projection and volume averaging operations 
to first transfer the forces experienced by the 
particles to the fluid and subsequently, extract 
the motion of the particle phase from the mo¬ 
tion of the fluid. With IBM, these operators 
are grid-dependent and are designed to pro¬ 
vide minimal resolution of the particles. The 
resulting flow due to the particles is typically 
computed using finite-difference or finite volume 
methods. With FCM, the operators are built 
around Gaussian functions which require more 
resolution than their IBM counterparts, but due 
to their smoothness, can be used in conjunc¬ 
tion with very accurate spectral solvers for the 
fluid flow. In addition, FCM aims to provide a 
higher resolution of particle finite size by includ¬ 
ing the operators that correspond to the second 
order effects of the rotlet and stresslet. The 
particle hydrodynamic interactions provide d by 
FCM have been studied extensivel} h^ l ^^ l ^^ and 
in many cases a successful comparison with SD 
results is found. 

To resolve Brownian motion with these meth¬ 
ods, the volume averaging operators can be used 
to obtain the random motion of the particles 
from the flow field generated by the fluctuat¬ 
ing stress. As the fluctuating stress is based 
on spatially uncorrelated white-noise, the ma¬ 
trix square root computation does not have to 
be performed. For the stochastic and fluctu- 
ating IBMs a nd fluctuating FCM, it has been 
showrPSESllIl explicitly that the resulting parti¬ 
cle velocity correlations satisfy the fluctuation- 
dissipation theorem. 

While the usage of fluctuating hydrodynam¬ 
ics has accelerated the computation of the ran¬ 
dom particle velocities, in the overdamped, or 
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Brownian dynamics limit where one can ignore 
both fluid and particle inertia, a seemingly- 
problematic Brownian drift term proportional 
to the divergence of the particle mobility matrix 
also needs to be accounted for even to achieve 
the correct equilibrium distribution of the sus¬ 
pension. 

Even though the drift term is identically zero 
for periodic suspensions when the point force 
solution is used to provide the mobility ma¬ 
trix, this is not the case when second-order 
corrections provided by the stresslets are in¬ 
cluded. This presented particular difficulty in 
developing fluctuating FCIVP^. While we were 
able to successfully and rapidly compute the 
correct random displacements with fluctuating 
FCM, the Brownian drift arising because of the 
stresslets limited any dynamic simulations we 
performed to only O(IO^) particles. In this 
work, we address this challenge directly and 
present a complete methodology for dynamic 
simulation with fluctuating FCM that accounts 
for the drift term due to the stresslets at a min¬ 
imal computational cost. 

Our approach relies on the construction of a 
time integration scheme t hat na turally accounts 
for the drift term. Fixmarp212fil showed that the 
effects of the drift term can be recovered with¬ 
out a direct computation by using a specific 
mid-point time integration scheme. This ap¬ 
proach, however, relies on the usage of random 
forces and torques, rather than random veloc¬ 
ities and angular velocities, making its imple¬ 
mentation with the stochastic and fluctuating 
IBMs and fluctuating FCM rather cumbersome 
and computationally expensiv^^. To overcome 
this difficulty, Delong et proposed an in¬ 
tegration scheme for the IBM using what they 
termed as random finite differencing (RFD). By 
randomly displacing and forcing the particles in 
a particular way, they showed that divergence 
of the mobility matrix due to point force hydro- 
dynamic interactions can be recovered without 
ever needing to perform the onerous computa¬ 
tions involved with Fixman’s method. 

In this study, we construct a time integra¬ 
tion for fluctuating hydrodynamics-based simu¬ 
lations of Brownian particles that include the 


second-order stresslet and rotlet corrections 
to the particle hydrodynamic interactions, for 
which RFD still incurs a significant computa¬ 
tional cost. We develop a midpoint time inte¬ 
gration scheme that we refer to as the drifter- 
corrector (DC) that requires the stresslet com¬ 
putation be performed only once per timestep, 
leading to the desired reduction in computa¬ 
tional cost. 

The DC works by first advecting the par¬ 
ticle positions a half timestep by the uncon¬ 
strained fluctuating flow field. At the midpoint, 
the complete calculation to determine the par- 
tice motion is performed, including the compu¬ 
tation of the stresslets. This time integration 
scheme with fluctuating FCM forms an efficient 
and complete method for the dynamic simula¬ 
tion of Brownian suspensions of moderate con¬ 
centration. In fact, we show that when the DC 
is used in conjunction with fluctuating FCM, a 
Brownian simulation requires just a single addi¬ 
tional unconstrained Stokes solve per timestep 
as compared to a deterministic FCM simulation 
of the same system. 

We provide an extensive validation of the 
DC by performing dynamic fluctuating FCM 
simulations of periodic suspensions and of a 
single particle between two slip surfaces. We 
show how to impose these boundary conditions 
with fluctuating hydrodynamics by modifying 
the spatial correlations of the fluctuating stress, 
rather than imposing it directly through the 
Stokes solver. This additional contribution pro¬ 
vides a framework to study Brownian parti¬ 
cles in thin films, allowing for direct compari¬ 
son with exper iment al results such as those on 
active particled^^EH addition, these simu¬ 
lations confirm our theoretical analysis of the 
scheme and shows that the DC is able to yield 
both the correct equilibrium distribution and 
the correct distribution dynamics as described 
by the Smoluchowski equation. 

We consider the canonical problem of a col¬ 
lapsing cluster of interacting colloidal parti¬ 
cles, allowing for comparison with results given 
by other methodologies and the demonstration 
that the higher-order stresslet corrections re¬ 
solved in fluctuating FCM yield quantitative 
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differences in the results. We show that, by 
avoiding additional stresslet iterations, the DC 
has a computational cost half that of forward 
RFD. 

Finally, using the DC with fluctuating FCM, 
we also perform simulations of suspensions of 
interacting Brownian particles to explore col¬ 
loidal gelation and percolated network forma¬ 
tion. We show that in addition to reproducing 
the results of other simulation techniques such 
as SD, fluctuating FCM with the DC allows for 
large-scale simulation of hydrodynamically in¬ 
teracting Brownian particles at a low computa¬ 
tional cost. 


II. EQUATIONS OF MOTION 


In this study, we will be considering a suspen¬ 
sion of Np spherical Brownian particles, each 
having radius a. The position of particle n is de¬ 
noted by Y”. Each particle may also be subject 
to non-hydrodynamic forces, F", and torques, 
r". In the overdamped or Brownian dynamics 
limit, w here b oth fluid and particle inertia are 
negligibl d^^ l ^^ l, the dynamics of the particles is 
described by the system of stochastic differen¬ 
tial equations 

dy = + M^^T) dt 

+kBTVy ■ M^^dt + dV (1) 

where y is the 3iVp x I vector that contains 
the position information for all particles, T is 
the 3iVp X I vector of forces on all of the parti¬ 
cles, and T is the similar vector for the torques. 
We note that because of their isotropic shape, 
there is no need to keep track of orientations 
for passive spherical particles and Eq. Q is 
sufficient to describe the dynamics of the sus¬ 
pension. Eor active or Janus spherical particles 
that possess an inherent swimming direction, 
one must also keep track of orientation and an 
additional equation similar to Eq. 0 for par¬ 
ticle orientation must also be considered. The 
vector dV is the incremental random velocity 
which, along with the incremental angular ve¬ 
locity, dVV, is related to the incremental QNp x I 


Wiener process, dB, through 


dV 

dW 


y/2kBTM^/'^dB 


( 2 ) 


where kB is Boltzmann’s constant and T is the 
temperature. Appearing also in Eqs. 0 and 
([^ is the 6Np x 6Np hydrodynamic mobility 
matrix 


M = 



Ajvr ■ 




( 3 ) 


as well as its submatrices and . The 
mobility matrix, which in general depends on 
the particle positions, provides the linear rela¬ 
tionship between the forces and torques on the 
particles and their resulting velocities and an¬ 
gular velocities. All the information about how 
the particles interact through the fluid is con¬ 
tained in the mobility matrix. The mobility ma¬ 
trix is determined by solving the Stokes equa¬ 
tions, 

Vp — TyV^u = 0 

V • u = 0 (4) 

that govern the flow induced in the surrounding 
fluid as the particles move through it. 


Eor methods such as Brownian dynamics or 
SD, the mobility matrix is constructed using 
the flow generated by a force multipole e xpan- 
sion in the Stokes equations and Faxen laws^IESl 
to extract the particle motion from the fluid 
velocity. In this paper, we adopt an alterna¬ 
tive, but paralle l app roach known as the force¬ 
coupling methocP^EU (FCM) that utilizes regu¬ 
larized, rather than singular, force distributions 
in the Stokes equations and replaces the Faxen 
laws by volume averaging operators. As noted 
by Delong el a/.^, this approach and the IBM 
share similar features and have distinct advan¬ 
tages over methods based on singular force dis¬ 
tributions, particularly when Brownian motion 
is included in the simulation. 

In the absence of Brownian motion, only the 
first term in Eq. Q, the mobility matrix mul¬ 
tiplied by the forces and torques, will be non¬ 
zero. The inclusion of Brownian motion gives 
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rise to the two additional terms - the incre¬ 
mental random velocity, dV, as well as an addi¬ 
tional drift known as Brownian drift. To sat¬ 
isfy the fluctuation-dissipation theorenP^, dV 
depends on the square root of the mobility ma¬ 
trix through Eq. ([^. This links the random 
motion of the particles with how they interact 
through the fluid. As we discuss in the next 
section, we can compute this term rather effi¬ 
ciently by volume averaging and constraining 
the fluid flow generated by a white-noise fluctu¬ 
ating stress. 

The Brownian drift term has a more subtle 
origin. It arises as a result of taking the over¬ 
damped limit of the Lan gevin dynamics where 
particle inertia is presenipEll. In order to suc¬ 
cessfully ignore the suspension dynamics dur¬ 
ing the short inertial relaxation timescale, one 
must include the Brownian drift term to obtain 
an SDE that is consistent with Smoluchowski’s 
equatiorPS. Note that since the particles con¬ 
sidered here are spherical, the Brownian drift 
term does not involve the divergence of the 
mobility matrix with respect to particle orien¬ 
tations, only particle positions. The purpose 
of this work is to introduce a new and study 
existin^islllllll time integration schemes that 
automatically account for Brownian drift. We 
pay particular attention to the case where the 
second-order correction to the particle mobility 
matrix is obtained by imposing a local rate-of- 
strain constraint on the flow. We show that for 
methods employing fluctuating hydrodynamics, 
the Brownian drift can be accounted for at the 
cost of a single additional unconstrained Stokes 
solve per time-step by using an appropriately 
designed time integration scheme. 


III. FLUCTUATING FCM 

To compute the terms in Eq. Q that de¬ 
pend on the mobility matrix, we utilize fluc¬ 
tuating ECM that is described and analyzed 
in detail in our pre vious worlP^. Fluctuating 
FCM combines FCIVP^I^ with fluctuating hy¬ 
drodynamics resulting in an efficient methodol¬ 
ogy to determine particle hydrodynamic inter¬ 


actions while simultaneously yielding the cor¬ 
rect random particle velocities that satisfy the 
fluctuation-dissipation theorem. In fluctuating 
FCM, the fluid velocity is given by the Stokes 
equations driven by a white-noise fluctuating 
stress, and a low-order, regularized multipole 
expansion representing the force the particles 
exert on the fluid. The motion of the particles 
is then found by taking volume averages of the 
fluid velocity. This process is similar to that 
of the stochastic and fluctuating IBMs, and as 
we demonstrate below, fluctuating FCM can be 
expressed in terms of projection, or spreading, 
operators and volume averaging, or interpola¬ 
tion, operators commonly used to describe IBM. 
We present fluctuating FCM using this frame¬ 
work to emphasize the connection between the 
methodologies, indicating that the time integra¬ 
tion schemes that we explore using fluctuating 
FCM in subsequent sections could also be used 
more widely. 

For fluctuating FCM, the fluid velocity is 
given by the following Stokes flow 

V-P-f 

+M^[T] + 1C^[S] 

V • u = 0, 

( 5 ) 

/C[u] = 0, (6) 

where the fluctuating stress, P, has the follow¬ 
ing statistics 

(^..(x)) = 0 (7) 

(Py (x)Pfc/(y)) = 2kBT{dik6ji + SiiSjk)S{x - y) 

( 8 ) 

with (•) denoting the ensemble average of a 
quantity. 

The non-hydrodynamic forces, P, and 
torques, T, on the particles, as well as the par¬ 
ticle stresslets, 5, the symmetric force-moment 
on each particle, are projected onto the fluid us¬ 
ing the linear operators A/”!, and which 
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are given by 

jt[j-] = ^F"A„(x) (9) 

n 

= ( 10 ) 

n 

/Ct[5] = ^ E S” • (V0"(^) + (V0„(x))^) . 

n 

( 11 ) 

Appearing in these expressions are the two 
Gaussian envelopes, or spreading functions, 

A„(x) = (27rai)-3/2e-l"-Y"l'/2-i (12) 

0„(x) = (27ra|)-3/2e-l’^-Y"l'/2-l (13) 

where cta and ere are related to the radius of 
the particles through cta = a/v^ and gq = 
a/ (d-v/rr) • Unlike the forces and torques 
which are typically set by external or inter¬ 
particle potentials, the stresslets arise as a result 
of the constraint on the flow given by Eq. 
and, consequently, need to be solved for as part 
of the general flow problem. The projection op¬ 
erators are the adjoints of the volume averaging 
operators J, N, and K. that are used to extract 
the particle velocities, V, angular velocities, W, 
and local rates-of-strain, £, from the fluid ve¬ 
locity and its derivatives, 



V = J[u] 

(14) 


yy = Af[u] 

(15) 


3 

I 

II 

(16) 

The expressions for these operators are most 
clearly expressed when they are restricted to 
particle n, 

V” = (J[u])„ = J uAnix)d^x 

(17) 

= (Af[u]) 

1 = ^ /u X V0„(x)d^x 

(18) 

E" = -(/C[u])„ 


II 

1 

tol 

u(V0„(x))’^ -k V0„(x)u^ 

d^x. 



(19) 


It is important to note that these definitions 
for the particle angular velocity and local rate 
of strain depart slightly from the standard FCM 
definitions provided by Lomholt & Maxej^H. As 
we show in Appendix using Eqs. and 

(19) ensures that, in general, the rate of work 


done by the particles on the fluid is correctly 
balanced by the viscous dissipation. These deh- 
nitions also ensure that the stresslets do no work 
on the fluid. If periodic or no-slip conditions 
are imposed on the bounding surfaces, integra¬ 
tion by parts of Eqs. and ( [l9| reveals that 
the particle angular velocities and local rates-of- 
strain reduce to the standard FCM definitions 
for these quantities which are the local volume 
averages of the fluid vorticity and rate-of-strain, 
respectively. With this in mind, the constraint 
Eq. ([^ insists the local rate-of-strain must be 
zero as is the case for a rigid particle, and the 
stresslets can be viewed as the Lagrange multi¬ 
pliers included to enforce this constraint. 

It is also important to note the dependence 


of the Gaussian envelopes, Eqs. (12) and (13) 


on particle position, Y". For a fixed flow field 
u, the values of V", rj", and E" are continuous 
functions of Y" and can be differentiated with 
respect to this quantity. We will see that this 
dependence is quite important and needs to be 
considered to design schemes that correctly up¬ 
date the particle positions when Brownian mo¬ 
tion is included in simulation. 

At this stage, it is useful to again note that 
the fluctuating and stochastic IBM share this 
framework with fluctuating FGM. The main dif¬ 
ferences between IBM and FGM are the choice 
of spreading function, and the only operators 
typically used with IBM are and its ad¬ 
joint. The higher-order correction to the hydro- 
dynamic interactions due to torques or stresslets 
resolved in FGM are not typically included with 
IBM. 


For FGM, the stresslets can be obtained us¬ 
ing the conjugate gradient procedure describe 
by Yeo and Maxey^, and once found, the so¬ 
lution to Eq. ([^ that satisfies Eq. (|^ is de¬ 
termined. For the simulations performed in 
the proceeding sections, approximately 10 con¬ 
jugate gradient iterations, each needing one 
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Stokes solve, are required to reach a residual 
of ||£|j = 7 • 10“®Z3o/a^, where 

||£||= max ||E"||2, (20) 

nG[l,iVp] 


volume averaging, and Stokes operators. For 
M'pcm’ expression is 

(23) 


and Dq = kBT/{&-KT]a) is the Stokes-Einstein 
diffusion coefficient for an isolated sphere. 

In designing, exploring, and constructing the 
time integration scheme that provides the Brow¬ 
nian drift term, we are particularly mindful 
of the computational cost associated with the 
stresslet computation, and in fact, seek to limit 
this computation to once per timestep. 


A. Mobility matrices and the fluctuation 
dissipation theorem 


Fluctuating FCM yields the deterministic 
and random velocities, corresponding to the 
first and third terms in Eqs Q without ever 
directly computing the mobility matrix. Al¬ 
though they are never computed explicitly, 
expressions for the mobility matrices can be 
determinecP^. Eor ECM, the mobility matri¬ 
ces provide the linear relationship between the 
deterministic particle velocities, angular veloc¬ 
ities, and local rates-of-strain and the forces, 
torques, and stresslets associated with the par¬ 
ticles. 


■ V ■ 


w 

£ 

— 


MV 


mV 


mV 


^FCM ■''''■FCM 

j^WS 


M 


'FCM 

ST 

FCM 


mV 


FCM 




M 
FCM 


T 
T 
S 
( 21 ) 

As shown in our previous worlP^, the submatri¬ 
ces can be expressed using the Gaussian spread¬ 
ing functions, Eqs. (12) and (13) and the 
Green’s function G(x,j^for the Stokes equa¬ 
tions. Eor example, the entries of the matrix 
MV;m linking particles n and m are given by 


^Fc’iV = J yA„(x)G(x,y)A„(y)(i^xdV- 

( 22 ) 

Similar expressions can be found for the other 
matrices in the appendix of our previous worlP^. 
As described by Delong et ali^, these matri¬ 
ces can also be expressed using the projection. 


where represents the inverse Stokes opera¬ 
tor. 

In the absence of the stresslets, the mobility 
matrix for ECM reduces to 


Mfcm = 


■'^'■FCM 

'^'FCM 


■''''FCM 

■''''FCM 


T 

T 


(24) 


If the stresslets are included, they can be found 
using Eq. (21). Since for rigid particles the local 
rates-of-strain are zero, we have £ = 0 and, in 
terms of the forces and torques, the stresslets 
will be given by 


S — —'R-Vcm i-^F^M^ + ■^Vcm'T') ■ (25) 


where we have written TVVcm — ■ 

From this expression for S, we find the stresslet- 
corrected mobility matrix is 


Mfcm-s 


f^VT f^VT 

■'^'fcm-s ■'^'fcm-s 

f^yvT f^wr 

■'^'fcm-s ■''''fcm-s 


(26) 


where 


\yiVT _ k/VT 

■I'''FCM-S — ■I'''FCM ~ ■''''FCM''-FCM-''''FCMj 

A/fVr _ f^vr f^vs -j^es ^sr 

FCM-S — ■''''fCM ~ ■''''fCM ''''FCM-''''fCM^ 

A/fWJG _ J^WT /^WS ^SS J^ST 

''''FCM-S — ''''fCM ~ ''''fCM '''-FCM''''fCM^ 

mWT _ Ji^wr /^ws ^ss J^sr 

FCM-S — ■''''FCM ~ ''''FCM '^FCM''''FCM^ 


(27) 


Again, while these matrices are never com¬ 
puted explicitly, we know that they exist. This 
allows us to draw a parallel between fluctu¬ 
ating FCM and traditional methods for sus¬ 
pended particles such as Brownian dynamic^^ 
and In addition, the expressions for 

the mobility matrices will prove useful in the 
analysis that demonstrates the time-integration 
scheme we propose in this study yields the cor¬ 
rect first and second moments of the incremen¬ 
tal change in the particle positions to first order 
in time. 

















The expressions for the mobility matrices can 
also be used to demonstrate that the random 
motion of the particles obtained using fluc¬ 
tuating FCM complies with the fluctuation- 
dissipation theoreirPS. If we consider the flow, 

Vp — = V • P 

V • u = 0 (28) 

/C[u] = 0, (29) 

the resulting particle velocities are then 

V = ^[u]. (30) 


By examining the velocity corrections explicitly, 
one can show^Sl that the particle velocities given 
by Eq. satisfy the fluctuation-dissipation 

theorenP^. For fluctuating FCM, we have that 

(VV^) = 2kBTM%M (31) 


when Eq. (29) is not enforced (i.e. 5 = 0) and 


(VV^) = 2kBTM%M-s (32) 


when the constraint Eq. ( |2^ is enforced. A 
similar result has been shown in th e con text of 
the stochastic and fluctuating lEM^SEH Thus, 
by simply applying the volume averaging opera¬ 
tors to the fluid flow induced by the white-noise 
fluctuating stress, we can obtain random parti¬ 
cle velocities consistent with the noise terms in 
the stochastic equations of motion. This is done 
without ever needing to explicitly form and de¬ 
compose the mobility matrix, leading to a great 
reduction in the computational effort needed to 
compute these terms. 


where one first solves the Stokes equations 
Vp'=-77V2u'= = (At)-5V-P'= 

V • u'= = 0, 

/C'=[u'=] =0, 


and then updates the particle positions using 

yk+i =3;fc + 

one does not account for ksT'Vy ■ 
and will not achieve the correct suspension dy¬ 
namics. Rather than resorting to computing 
the drift term explicitly, however, its effects 
can be successfully incorporated into a simu¬ 
lation by using an approp riately designed time 
integration schem#2I2212Il_ jn this section, we 
discuss these schemes, and building from the 
ideas used in their construction, we introduce 
a new scheme called the drifter-corrector (DC). 
The DC has the particular advantage that when 
constraints are imposed on the flow to recover 
stresslet corrections for the particle hydrody¬ 
namic interactions, they need only to be im¬ 
posed once per timestep. This leads to a non- 
negligible reduction in computational cost with 
respect to other schemes. Eor the DC, the to¬ 
tal additional computational cost per timestep 
for fluctuating FCM with respect to its deter¬ 
ministic counterpart is a single Stokes solve per 
timestep. 


A. Fixman’s method 


IV. TIME INTEGRATION 

While fluctuating FCM provides the deter¬ 
ministic and random particle velocities corre¬ 
sponding to the first and last terms in the equa¬ 
tions of motion Eq. ([^, the Brownian drift, 
fcs TVy ■ would also need to be ac¬ 

counted for to advance the particle positions in 
time. Therefore, simply W applying the Euler- 
Maruyama (EM) schem^^ to fluctuating FCM 


In the late 1970’s, FixmaiP^*2fil developed a 
midpoint integration scheme to account for the 
drift term. The key to recovering the drift 
is that the scheme employs the same random 
forces at time levels tk and t/c+i /2 while using 
updated values for the particle positions at level 
tk+i/ 2 - Applying Fixman’s method to fluctuat¬ 
ing FCM yields the following scheme: 

1. Compute the random velocities, and 
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angular velocities, 

Vf - 7?V^u'= = (At)-5 V • 

V • = 0, 

K.’^lu^] =0, 

= J[u^] 

= AA[u'=] 


2. Compute random forces, , and torques, 
r'^, from and 


'F>^' 

A /f ~1 

■ yk ■ 

'j-'k 

= -^FCMS 

Wk _ 


(33) 


3. Solve the predictor constrained Stokes 
problem 

Vp'^ - TjV^u'^ = + F^\ 

V • = 0 

lC^[u^] = 0 

(34) 


The subscripts k and fc + 1/2 indicate whether 
the operators are evaluated at the positions 3^^, 
or 3 ;^+i/ 2_ This scheme provides the first and 
second moments of the increment up to first or¬ 
der in At. We note that the approximation of 
the deterministic aspects of the motion can be 
improved to second order by using and 

q-k+i /2 instead of F^ and in Step 5. 

While Fixman’s scheme does avoid a direct 
calculation of the Brownian drift term, it does 
introduce the need to solve a large linear system 
in Step 2 to find the random forces and torques. 
This system can be solved using conjugate gra¬ 
dient methods, with each iteration requiring the 
Stokes equations to be solved. Simulations with 
N = 183 particles corresponding to a volume 
fraction of ^ = 0.1 required approximately 25 
iterations to find the random forces and torques 
with the residual reduced to le-4 times its origi¬ 
nal valu^^. This computation, along with hav¬ 
ing to compute the stresslets twice per timestep, 
severely limit the scale at which Brownian simu¬ 
lations could be performed using the fluctuating 
FCM. 


B. Random Finite Differencing 


4. Move particles to the midpoint, 

yk + l/2 ^yk ( 35 ) 

5. Solve the constrained Stokes problem at 
the midpoint, 

^pk+l/2 _ ^y2^fe-ei/2 ^ jP,k+l/2yj:k 

^^t;/c+l/2[5/c+l/2] 

V • u'=+i/2 = 0 
/C'=+l/2[u'=+l/2] = 0 

(36) 


6 . Update the particle positions 

yk+i ^yk 

( 37 ) 


To help overcome this challenge, Delong et 
introduced the random finite difference 
(RFD). RFD takes advantage of the fact that 
the random forces and torques and the result¬ 
ing displacements used in Fixman’s method are 
in fact one of many possible choices that could 
be used to account for Brownian drift. Specif¬ 
ically, they show that using random displace¬ 
ments 5Ajy, as well as randomly particle forcing 
AF/6, the divergence of a mobility matrix, A4, 
can be approximated from 

^{M{y + 5Ay) AF - M{y)AF) = VyM 

+o(S) 

(38) 

provided that (AJ^AJ") = X. Therefore, one 
can simply choose AJ^ = AF = where ^ is 
a 3JVp X 1 vector of independent Gaussian ran¬ 
dom variables with zero mean and unit variance. 
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This eliminates the need to solve any linear sys¬ 
tem. The concept of random finite differencing 
can also be extended to higher-order accuracy. 
For example, the central random finite differ¬ 
ence is given by 

^{M{y + -M(y- 51/2)0 

= VyM+0(S^) 

(39) 

Thus, in general, a weakly first-order accurate 
scheme that accounts for the Brownian drift 
term can be constructed by adding an RFD to 
the Euler-Maruyama scheme. Applying this to 
fluctuating FCM, we have: 

1. Solve the constrained Stokes problem 

V/ - = (At)-5 V • 

V • = 0 

/C'=[u'=] = 0 (40) 

2. Solve the constrained Stokes problem with 
disturbed particle positions but without 
the fluctuating stress tensor, 

V/+‘5 - 

V • = 0 

/C'=+'5[u'=+'^] = 0 (41) 

3. Add the corresponding particle velocities 

and move particles to the next timestep, 

y>‘+^ = y^ + At 

+Atj'^+^[u'^+^] (42) 

where the superscript k + S indicates that the 

operators are evaluated at y'^ + 6$. We em¬ 
phasize that for passive spherical particles, we 
only need the divergence of the velocity-force 
mobility matrix. We, therefore, only need to 
consider random displacements and forces with 


the RFD. A similar scheme can be constructed 
using central RFD to account for the Brown¬ 
ian drift. Comparing the resulting scheme with 
Fixman’s method, we immediately see that us¬ 
ing RFD eliminates the need to compute any 
random forces or torques, making this approach 
much more suited for fluctuating FCM. 


C. Drifter-Corrector 

While the usage of RFD provides a clear ad¬ 
vantage over Fixman’s method for fluctuating 
hydrodynamics-based methods, we see that it 
does require that the stresslets be determined 
twice per timestep (three times per timestep if 
using central RFD). The cost of performing a 
Brownian simulation is then at least double that 
of a deterministic simulation of the same sys¬ 
tem. For fluctuating FCM where the conjugate 
gradient method is used to find the stresslets, 
the additional cost of including Brownian mo¬ 
tion would then be approximately 10 Stokes 
solves per timestep. 

In an effort to mitigate this computational 
cost as much as possi ble, w e build on the ideas 
introduced by FixmarP^^ and Delong et 
and construct a mid-point scheme that accounts 
for Brownian drift at the cost of one Stokes solve 
per timestep. We will refer to this scheme as the 
drifter-corrector (DC). Specifically, the DC is: 

1. Solve the unconstrained Stokes problem 

V/ - = (At)-5 V • 

V • u*’ = 0 (43) 


2. Move particles to the midpoint, 


yk+i/2 ( 44 ) 
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3. Solve the constrained Stokes problem, 
V/+1/2 _ ^^2^k+l/2 ^ (At)-5 V • P'^ 

j-t;fe+i/2[j-fc+i/2j 

_l__;yt;fc+i/2j7-fc+i/2] 

^^t;fc+i/2[5fc+i/2] 

V • U*^+l/2 = 0 

/C''+^/2[u'=+i/2] = 0 (45) 

4. Update the particle positions 

yk + l ^ yk 

+Atil+v'^) 

xj^fe+i/2[ufe+i/2], (46) 


The scalar factor, is related to the diver¬ 
gence of the unconstrained particle velocities, 
\5t = (^'=[u'=])„, through 

= (47) 

n 


For fluctuating FCM, we compute u* directly 
by evaluating 


V 


k 



u^A„(x)d^x(48) 


but it can also be computed using finite- 
differencing or even RFD. 

We show explicitly in Appendix that the 
expansions of the DC for the first and second 
moments of the increment 
are 

{Ay^} = AtMpQj^f_gF^ + 

+AtkBTVy ■ MYcii-s + 0(At2), 

(49) 

and 

(A3;'=(A3;'=)^) = 2kBTAtMYcM-s 

+0{AY, (50) 


The DC, therefore, first moves the particles 
a half timestep using the particle velocities ob¬ 
tained from the unconstrained fluctuating flow 
field. Then, using the updated positions, but 
the same realization of the fluctuating stress, 
the particle forces and torques are projected 
onto the fluid, and the local rate-of-strain con¬ 
straint is imposed. The particle positions are 
then updated using the resulting particle veloc¬ 
ities. Thus, for the DC, the stresslets need to be 
computed only once per timestep, and the only 
additional cost is the single Stokes solve to de¬ 
termine the unconstrained fluctuating flow field. 
A second-order approximation of the determin¬ 
istic terms may be achieved by also including 
particle forcing (torques, forces, and stresslets) 
in Eq. (43), but this would come at the cost of 


an additional stresslet iteration per timestep. 

We note that similar ideas were employed by 
Delong et to construct a midpoint scheme 
using RFD. In their construction, however, they 
relied on a specific decomposition of the mo¬ 
bility matrix that is not applicable in the case 
where local rate-of-strain constraints are en¬ 
forced. Nevertheless, our analysis reveals that 
this similar technique extends to the case where 
the stresslets are present. 

If we have periodic boundary conditions, or 
if u • n = 0 pointwise on the boundary, where 
n is the unit normal to the boundary, an anal¬ 
ysis of the scheme using continuous spatial op¬ 
erators reveals that the correct first and second 
moments can be achieved with = 0. This 
condition holds for no-slip or no shear on the 
fluid boundaries. As shown in Appendix [B} this 
is derived from the condition that V y" • U„ = 0 
for all particles. In this case, Eq. ( [4^ simplifies 
to become 

yk+l ^yk ( 54 ) 


It should be emphasized that these arguments 
based on continuous spatial operators may not 
hold in practice when the Stokes equations and 
the spreading and volume averaging operators 
are discretized. We may not, therefore, nec¬ 
essarily set = 0 without incurring an addi¬ 
tional small error due to spatial discretization. 
Delong et have shown for the fluctuating 


respectively. 
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IBM that even in the case of periodic boundary 
conditions, spatial discretization does give rise 
to this small error, but it can be eliminated by 
including an RFD of the spreading operator as 
part of the time integration. For the DC, this 
would correspond to including the term. In 
Appendix [B] and Section [V B[ we explore these 
errors for fluctuating FCM. We hnd that for 
periodic boundary conditions, the spatial dis¬ 
cretization does not introduce any additional 
errors for fluctuating FCM and we may use the 
scheme with = 0. We do find, however, that 
for no-shear conditions on the boundaries, the 
numerical integration of the Gaussian functions 
introduces this small, but non-negligible, error. 
Thus, we retain the terms in our simulations 
where we apply no-shear boundary conditions. 


V. NUMERICAL STUDIES 

To demonstrate the performance of the DC 
with fluctuating FCM, we perform simulations 
of particulate suspensions under both dynamic 
and equilibrium conditions. These simulations 
confirm the results of our theoretical analysis of 
the DC, and show, in practice, that it is able to 
produce the correct microstructure, dynamics 
and final equilibrium states for distributions of 
particles. Our results also show that with the 
DC, fluctuating FCM can be used for large-scale 
simulation of Brownian suspensions even when 
higher order corrections to the hydrodynamic 
interactions are included. 


A. Spatial discretization 


In our simulations, we mainly consider pe¬ 
riodic boundary conditions and use a Fourier 
spectral method with fast Fourier transforms to 
solve the Stokes equations, taking advantage of 
the highly scalable MPI library PSDFFTI^. We 
give here a summary of the main steps of the 
discretization scheme as a detailed description 
is provided in our previous worlP^. 

For the case where each side of the domain 
has length L, we use M grid points in each 


direction. The total number of grid points is 
then Ng = and the grid spacing is given 
by Aa; = L/M. The position of a grid point 
in a given direction is then Xa = aAx for 
a = 0 ,..., M — 1. At each grid point, the en¬ 
tries of the fluctuating stress are independent 
Gaussian random variables. Since the fluctuat¬ 
ing stress P is symmetric, at each grid point 
we generate six random numbers based on a 
Gaussian distribution with zero mean and unit 
variance. We then multiply the off-diagonal en¬ 
tries by y/2kBT/{Ax^) and the diagonal ones 
by 2yJkBT/{Ax^). 

To the fluctuating stress, we add the FGM 
force distribution, 

fp’CM(x) = 

+AfnT](x)+/C^[5](x) (52) 


evaluated at the grid points, x = 
[xa xp x^Y', taking also A„(x) = 0 and 
0„(x) = 0 for |x — Y”! > 3a. To ensure suffi¬ 
cient resolution of ipCM, we set ae/Ax = 1.5 
and a^/Ax = 1.86. After taking the discrete 
Fourier transform (DFT) of the total force 
distribution, we compute the fluid velocity in 
Fourier space 


U (/t Q, , , /t .y ) 


1 _ kk 


r^jkp \ |kp^ 

X (ik- {At)-^/^'P{ka,k0,kY 

+fFCM(^Q, fc/3, Y) 


(53) 

where k = [ka kp k^ ]^. We set u(k) = 0 if 
|k| = 0 or if k • = Mtt/L for z = 1, 2, or 3. 

After taking the inverse DFT to obtain the fluid 
velocity at the grid points, the velocity, angular 
velocity, and local rate-of-strain for each parti¬ 
cle is computed by applying the spectrally accu¬ 


rate trapezoidal rule to equations Eqs. (17 - 19), 
where we again set A„(x) = 0 and = 0 

for |x-Y"| > 3a. 

For the case where the stresslets are ignored, 
we have 5 = 0 and the Stokes equations need 
only to be solved once per time-step to obtain 
the particle velocities. If the stresslets are in¬ 
cluded, however, they must be solved for, and to 
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do this, we employ the iterative conjugate gra¬ 
dient approach described in Yeo and Maxey^. 


B. Brownian particles in a periodic domain 

In this section, we perform two simulations to 
test the DC with fluctuating FCM in periodic 
domains. We first simulate the free diffusion of 
non-interacting particles to quantify the effect 
of spatial discretization of the operators. This 
same test was performed by Delong et for 
IBM. We then simulate a hard-sphere suspen¬ 
sion to demonstrate that the DC successfully ac¬ 
counts for the Brownian drift term that arises 
due to the stresslets. We show that with the 
DC, the correct radial distribution function is 
reproduced, while neglecting the Brownian drift 
term by using the EM scheme does not provide 
the correct suspension microstructure. We note 
here that as we have periodic boundary condi¬ 
tions, we take = 0 in our DC simulations. 

1. Free difFusion 

For IBM and FCM without stresslets, the 
translational invariance associated periodic 
boundary conditions indicates that = 

0. Delong et however, showed that while 

this holds true for the spatially continuous for¬ 
mulations of these methods, the spatially dis¬ 
crete operators used in simulations may not pre¬ 
serve translational invariance, and hence the di¬ 
vergence of the mobility matrix may be non¬ 
zero. Indeed for the IBIVP^ , spatial discretiza¬ 
tion leads to a non-zero divergence and when 
EM is used to integrate the equation of motion, 
variations in the particle distribution are found 
at the level of the grid spacing even though the 
distribution should be uniform. 

To quantify this effect for FCM, we let 522 
non-interacting particles diffuse in a 3D periodic 
system for t/toa = 1^0 Brownian times where 
/Dq is the characteristic time for a par¬ 
ticle to diffuse one radius and Dq is the Stokes- 
Einstein diffusion coefficient. The stresslets are 
not included in these simulations, making this 


test case identical to that performed by Delong 
et al^. We perform 20 independent simula¬ 
tions, each run using the very small time-step 
At = A.\e-A{Ax^/ Dq) where Ax is the grid 
spacing. This small value of At is chosen to 
ensure that any systematic variations in the dis¬ 
tribution are due to the spatial discretization of 
the operators. Fig. shows the depth-averaged 
normalized particle distribution over a single 
grid cell for FCM simulations using both EM 
and the DC. To allow for a quantitative com¬ 
parison with Delong et we have set the 

limits of the color axis exactly as it is in Fig. 3 
of their work. Apart from infinitesimal statis¬ 
tical fluctuations, there is no clear bias in the 
distributions for either integration scheme. This 
means that our spatial discretization preserves 
translational invariance within a tight tolerance. 
This property is likely due to the smoothness of 
the Gaussian kernels used with FCM, as well 
as the spectral accuracy of the Fourier method 
used to solve the Stokes equations. 


2. Hard-sphere suspension simulations 

When stresslet corrections to the particle hy¬ 
drodynamics interactions are included, the di¬ 
vergence of the mobility matrix is no longer 
zero for periodic systems. We show here that 
time integration with the DC recovers the drift 
term due to the stresslet corrections in FCM 
and show that neglecting the drift term by us¬ 
ing EM leads to an erroneous characterization 
of the suspension microstructure. 

To do this, we simulate a suspension of 174 
particles in a periodic domain that interact via 
a hard-sphere potential. The hard-sphere po¬ 
tential is included by rejecting timesteps where 
the particles are observed to overlap. The vol¬ 
ume fraction is taken to be </) = 0.1. We per¬ 
form these simulations using both EM and the 
DC to advance the particle positions in time. 
For each time integration scheme, we perform 
20 independent simulations, taking each simu¬ 
lation to a final time of t/tjj^ = 125. Fig. 
shows the radial distribution function g[r) ob¬ 
tained by averaging over the EM simulations 
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Figure 1. Normalized depth-averaged equilibrium distribution in a grid cell for simulations using (a) EM, 
(b) the DC 


and the DC simulations, as well the values given 
by the Percus-Yevick approximatiorP^. The DC 
simulations are in very good agreement with 
the Percus-Yevick approximation for all values 
of r, showing that the drift term is accounted 
for in these simulations. The EM results, how¬ 
ever, deviate from the Percus-Yevick approxi¬ 
mation for r < 2.6a and clearly overestimate 
g{r) at contact. The separations r < 2.6a cor¬ 
respond to the range over which the stresslets 
strongly modify the mobility matrix. These 
results clearly indicate that even for periodic 
boundary conditions, stresslet corrections to the 
hydrodynamic mobility yield a non-negligible 
drift term, which is correctly accounted for by 
the DC. 


C. Brownian particles between two slip surfaces 

We examine the Brownian dynamics of a 
spherical particle between two slip surfaces at 
z = 0 and z = L/2. On these surfaces, the 
boundary conditions for the flow are given by 
u • z = 0 and (I — zz"^) Vu = 0. The divergence 
of the mobility matrix is non-zero for these con¬ 


ditions and the drift term must be properly ac¬ 
counted for in simulation. Even though we im¬ 
pose periodicity or have u • n = 0 on the bound¬ 
aries, it is important to include in the time 
integration to eliminate errors that arise when 
the Gaussian functions overlap the boundaries 
(see Appendix [b| . 


1. Image system and fluctuating stress for slip 
surfaces 

Rather than enforcing these boundary con¬ 
ditions directly through the numerical solver, 
we can create the slip surfaces by introducing 
the appropriate image system for both the FCM 
particle force distribution, and the fluc¬ 

tuating stress, P. The flow due to the force 
distribution and its image can then be found 
using the Fourier spectral method described in 
Section IV A1 

In our simulations, the fluid domain is given 
by D = [0, L]^ X [0,^2]. We impose the slip 
conditions at z = 0 and z = and periodic 
boundary conditions on the other boundaries 
such that u(0, j/, z) = u(L, y, z) and u(x, 0, z) = 
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u(x, L, z). The force distribution due to particle 
forces is given by Eq. ( [^ . To enforce the slip 
conditions at z = 0 and z = L^, we first extend 
the fluid domain to ft* = [0,L]^ x [0, 2L^] and 
utilize the modified FCM force distribution and 
its image, 


Figure 2. Radial distribution function for a semi¬ 
dilute suspension (</> = 0.1) of hard-spheres in a 
periodic box. The blue line with square markers 
corresponds to simulations with EM time integra¬ 
tion, the black line with circles corresponds to the 
DC, and the red line is given by the Percus-Yevick 
approximat iorP^. 


fFCM(x) 


^t[j-](x) + Qt['^(x)+/Ct[5](x) , XG [0,L]2 X [0,L,] 

0 , xG [0,L]2 x]L„2L,[ 


(54) 


fFCM,im/\ _ j 0 , X G [0, L]^x]0,Lz[ 

W \ (J - 2zz^) f^^^(X) , X G [0, L]2 X [L„ 2L,] ' 


(55) 


where X = x — 2(x-z — Lj)^. The total forcing 
term that will now appear in the right hand 
side of Eq. (|53[ ) is given by the summation 
fFCM _|_ fFCM,im^ jg determined 

over the domain 11* with the periodic boundary 
conditions in all directions. The combined ef¬ 
fects of the image system and periodicity in the 
z-direction yields the desired slip conditions at 
z = 0 and z = L^- Once the fluid flow is deter¬ 
mined, the particle velocities, angular velocities, 
and local rates-of-strain are computed using the 
resulting flow restricted to the domain O. 

In introducing the slip surfaces, one must 


r 


make a choice as how to modify the Gaussian 
functions when the particles get close to the 
boundaries. In this work, we simply truncate 
any part of the Gaussian functions that ex¬ 
tends beyond the domain O. This is already 
reflected in our definition of the force distribu¬ 


tion, Eq. (54) and its image, Eq. (55). We use 


the same truncated Gaussians to compute the 
particle velocities, angular velocities, and local 
rates-of-strain, thereby preserving the adjoint 
properties of the operators. For particles in con¬ 
tact with the boundary, the truncated volume 
in Eq. 0 is approximately 3.8% of the total 
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and for Eq. (191 it is only 1.4%. We note, that 


more sophisticated ways to modify the Gaus¬ 
sian envelopes have been explored and are im¬ 
plemented in FCM elsewher^^. 

Along with the FCM particle force dis¬ 
tributions, the fluctuating stress must also 
have the appropriate symmetries to satisfy the 
fluctuation-dissipation theorem when the slip 

boundaries are present. In Appendix we 

show explicitly that this is achieved by having 

(P.,(x)) = 0 (56) 

(Pij (x)Pfe/(y)) = 2kBT- y) 

+2kBTV,,ki5{^ - Y) (57) 

for X, y G n* where 

^ijkl — ^ik^jl “t" ^il^jk (^^) 

^ijkl — 4“ (^^) 

with = 5jk - 2 d 3 jS 3 k, and 


Y = y - 2{y ■ z - 


(60) 


The symmetrized fluctuating stress Eq. (56 )- 
( [^ is discretized at each grid point with in¬ 
dependent Gaussian random variables whose 
statistics are 


(Tlij (^Q , ^/3, ^ 7 )) — 0 


(Tltj , 3^/3 , ^ 7 ) 7 ^/ (^a , ^/3 j ^K,)) — 
(Ax)3 ^ijklSjK + (Ax)3 


(61) 


(62) 


stresslets to obtain ^{z) and 

respectively. The values of 2kbT^^^~^ 


We determine the coefficient n±{z) for FCM 
by applying a unit force F = i on an iso¬ 
lated sphere and measuring the resulting ve¬ 
locity Vz(z) at 400 equi-spaced values of z be¬ 
tween z = a and z = — a. We perform these 

simulation both with and without the particle 

- . ^ n 

and 2kbT^^^^ (z) are provided in Fig. We 
see that both with and without the stresslets, 
the value of this mobility coefficient decreases as 
the particle approaches the slip surface. The ad¬ 
dition of the stresslet results in the mobility co¬ 
efficient depending more strongly on z and fur¬ 
ther reduces the mobility coefficient by approx¬ 
imately 30% near the boundaries. We have per¬ 
formed similar computations for (z) by tak¬ 
ing F = i, but have not included this data for 
brevity. We also note that given the symmetries 
associated with a single slip surface, the exact 
Stokes flow for a rigid sphere near a slip sur¬ 
face is equivalent to an appropriate two-sphere 
problem. A comparison of FCM with these two- 
sphere problems has been addressed by Lomholt 
and MaxejS^, showing that FCM provides ac¬ 
curate results for a wide range of separations. 

To ensure that the fluctuation-dissipation 
theorem for the particles is satisfied when the 
statistics for the fluctuating stress are given by 
Eq. (56l-(57), we compute the component of 


the random particle velocity Vj_(z) normal to 
the slip surface. We then can check that its 
correlations satisfy 


where the index k*™ = mod (M — k,M). We 
remind the reader that the indices a, /3,7 , k go 
from 0 to M — 1, where M is the number of 
discretization points in a given direction. 

2. Single particle mobility 

With the slip boundaries present, the mobil¬ 
ity matrix for a single particle has the form 

=/rj_(z)zz^-I-^11 (z)(I — zz^) (63) 

where the mobility coefficients p.||(z) and /i_L(z) 
depend on the distance from the slip surfaces. 


To efficiently compile a large number of sam¬ 
ples to accurately determine {z), we di¬ 

vide the domain in the z-direction into 200 eq- 
uispaced parallel planes. Over each plane, we 
distribute 50 particles that do not interact with 
one another. For the stresslet case, particular 
care is taken to ensure that any interactions 
are removed. This is done by utilizing directly 
the stresslet-corrected mobility coefficients mea¬ 
sured from simulations of an isolated particle. 
The procedure is detailed in Appendix 
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Figure 3. A comparison between the wall-normal 
mobility coefficient and the autocorrelation of the 
wall-normal particle velocity. 



z = 2Lj = 0 

'y = M = 0 


z = Lz 
7 = M/2 


2 = 0 

7 = 0 


Figure 4. Sketch showing a particle (and its image) 
in the slip channel. The black dashed line repre¬ 
sents the periodic boundaries conditions. The in¬ 
dex 7 runs over the z-coordinate of the grid. The 
potential [/(z) provides the interaction between the 
particle with the slip surfaces and is given by Eq. 
(651. The force on the particle is Fu{z) = —dU/dz. 


For 10000 realizations of the fluctuating 
stress, we compute V±{z) for each of the par¬ 
ticles. For each fixed value of z, we compute 
(z) by averaging over the realizations and 

the 50 particles at that particular value of 0. 
Fig. [3] shows the correlations of the particle 

normal velocity, At/Vl)(z), and the normal 


mobility coefficient in the channel, 2kBT^±{z). 
We see that for both cases, with and without 
the stresslets, the fluctuation-dissipation rela¬ 
tion for the particles, Eq. (641, is satisfied. 
Though not shown, we have also performed the 
same check for the parallel random velocities 
and a similar strong agreement was found. 


3. Equilibrium distribution between two slip 
surfaces 


To demonstrate their RFD approach, Delong 
et performed simulations of one and 100 
spherical particles between two no-slip surfaces 


at z = 0 and z = and subject to the potential 


U{z) = 


yz-Rreff 


, Z Rref ) 


, Z > Lz Rref j 


. 0 , otherwise. 

(65) 

When the particle comes within a distance Rref 
from the wall, it experiences linear, repulsive 
force whose strength is governed hy k. If the 
time integration scheme correctly accounts for 
the Brownian drift term, the equilibrium dis¬ 
tribution of the particles will be given by the 
Boltzmann distribution. 


P{z) = Z ^ exp 


uyy 

kBT_ 


( 66 ) 


where Z = 


exp 


Ujz) 

kBT 


dz. 


We consider a similar problem here, now sim¬ 
ulating the dynamics of a spherical particle be¬ 
tween two slip surfaces, but still subject to the 
potential U{z). A sketch of this simulation is 
provided in Fig. Even though we have slip 
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boundary conditions on the walls, the distribu¬ 
tion of particles should still be described by Eq. 
([66 ). We perform these simulations both with 
and without the stresslets and using both cen¬ 
tral RED and the DC. We also integrate the 
equations using the Euler-Maruyama scheme, 
Eq. ( |33| , which should instead yield the biased 
distributiorfi^, 


Pb{z) = Z^^exp 


U{z) -I- /cBTln(^_L(z)) 
keT^ 


(67) 

In our simulations, we have set = 9.7a, Lx = 
Ly = 2Lz, Rref = l-4a and k = 24kBT/(Ax)"^. 
For the RED simulations we set S = 10“®Aa::. 
To obtain the equilibrium distribution, each 
simulations is performed with Np = 500 non¬ 
interacting particles and run to time t/to = 

a2 

391.60, where to = is the time re- 

kBT 

quired for the particle to diffuse one radius in 
the z-direction based on the z-averaged mobil¬ 
ity coefficient 


M-L 


rLz —a 


Lz — 2a , 


/j,±(z)dz. 


( 68 ) 


To determine P{z) from the simulations, a 
histogram of the particle positions in the z- 
direction is compiled during the time window 
t/to^ = 6.526 — 391.60. The results from the 
simulations without the stresslets are shown in 
Fig. |5a] wh ile those with the stresslets are given 
in Fig. |5b| In both figures, we see that central 
RED and the DC recover the correct Boltzmann 
distribution, Eq. (66). For completeness, we 
have also run the DC simulations with = 0. 
We see that since the potential prevents the par¬ 
ticles from getting close to the boundaries, the 
= 0 simulations still provide the correct dis¬ 
tribution. The Euler-Maruyama scheme, how¬ 
ever, yields the biased distribution, Eq. (67), 
as the Brownian drift, fc^Td/ij_/dz, is non-zero. 
We see also that compared to the stresslet-free 
case, when the stresslets are included the biased 
distribution exhibits higher peaks close to the 
boundaries. This is a direct result of the greater 
reduction in mobility near the boundaries when 
the stresslets are included, see Fig. One can 


then expect that using the correct integration 
scheme becomes even more important in sim¬ 
ulations where the singular lubrication interac¬ 
tions with the boundaries are resolved. 


4. Distribution dynamics 


Along with characterizing the equilibrium 
distribution, we perform simulations to confirm 
that fluctuating FCM with the DC also yields 
the correct distribution dynamics. For non¬ 
interacting particles, the dynamics of the distri¬ 
bution is described by the Smoluchowski equa¬ 
tion 


dP 

iH 


d_ 

dz 


^ii_FuP - kBT^Mx{z) 


dP' 

dz 


(69) 


for the distribution P{z,t) subject to no flux 
conditions 


H±FuP - kBTfj.± 


dP' 

dz 


= 0 

z—Q,Lz 


(70) 


at the slip surfaces (z = 0, z = Lz). The ad- 
vective flux term is proportional to the force 
Fu = —dU/dz associated with the potential U. 
We solve Eq. ( |6^ using a first-order finite vol¬ 
ume solver where the advective fluxes are calcu¬ 
lated using an upwinding scheme and the diffu¬ 
sive terms with a central scheme. We have val¬ 
idated the solver against analytical solutions of 
the heat and transport equations. The mobility 
coefficient fj,± (z) is determined by interpolating 
the values from our FCM simulations, see Fig. 

to the finite volume grid points. 

We obtain the distribution dynamics from 
stresslet-corrected fluctuating FCM simulations 
by averaging the histogram time dynamics over 
20 independent simulations each with 500 non¬ 
interacting particles. For these simulations, the 
initial particle positions are generated by dis¬ 
tributing their positions uniformly in the center 
of the domain in the region z G [3Lz/8-, 5L2/8]. 
We compare these results with our numerical 
solution of Eq. (69) using the distribution ob¬ 
tained from the fluctuating FCM simulations at 
t = 0 for the initial condition. 
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Figure 5. The equilibrium distribution for a Brownian particle subject to the potential U{z). The distri¬ 
bution from simulations is obtained by averaging 10 simulations of 500 non-interacting particles over the 
time interval t/toa ~ 6.526 — 391.60. —©- : DC without : DC with v'^, —$— : central RFD 

with 5 = 10“®Aa:, —#— : Gibbs-Boltzmann distribution Eq. (66l, — B —: Euler-Maruyama Scheme, —■ —: 
Biased Gibbs-Boltzmann ditribution Eq. (|67[). (a) Without stresslets, (b) With stresslet corrections. 


Fig. [6] shows the dynamics of the distribu¬ 
tion obtained from the fluctuating FCM simu¬ 
lations with the DC time integration, as well 
as the finite volume solution of Eq. (691. We 
see that the distributions given by both meth¬ 
ods match as they evolve from the initial to the 
equilibrium state reached at time t/to^ = 6.4. 
This confirms further that the DC recovers the 
dynamics described by the Smoluchowski equa¬ 
tion, Eq. (69), as desired. 


D. Colloidal gelation and percolation 

An important application of the methodology 
presented in this paper is the large-scale simula¬ 
tion of interacting colloidal particles. Depend¬ 
ing on the interactions between these particles, 
colloidal suspensions c an ta ke on both solid- 
and fluid-like propertie^^EZl When the parti¬ 
cles interact via an attractive potential, the sus¬ 
pension can transition from a fluid-like state to 
that of a solid through the process of gelatiorP^l. 


This transition depends not only on the poten¬ 
tial, but also on the volume fraction occupied 
by the particles. Recent numerical studies have 
shown that the dynamics of this aggregation 
process will also depend on the hydrodyna mic 
interactions between the particleJi^l2lIllllZI_ jn 
this section, we use fluctuating FCM with the 
DC to simulate of the aggregation and gela¬ 
tion of interacting particles, focusing particu¬ 
larly on the role of hydrodynamic interactions 
on the resulting structures. We co mpare our re¬ 
sults with previous studie^l^ ^^ * ^^ l, demonstrat¬ 
ing that fluctuating FCM with the DC provides 
an effective and computationally efficient ap¬ 
proach for studying Brownian suspensions. 


1. Collapsing icosahedron 

A simple example that highlights the role of 
hydrodynamic interactions on colloidal aggrega¬ 
tion is the collapse of a small cluster of Brow¬ 
nian spherical particles as originally studied by 
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Figure 6. The evolution of the particle distribution 
during the time interval tjtoa = 0 — 6.8037. The 
results from the particle simulations are achieved by 
averaging over 20 simulations, each with 500 non¬ 
interacting particles. The stresslet corrections are 

included in these simulations.-: fluctuating 

FCM with the DC, -: solution to the Smolu- 

chowski equation (691 - (701. 


Furukawa and TanakaP^. Specifically, they con¬ 
sidered Np = 13 Brownian particles initially 
at the vertices of a regular icosahedron with 
edge length S.OSAx. The particles interact via 
a modified Asakura-Oosawa potential 


fc — Hj) r ,r < D 
UAo{r) = 'I c (T)| - l/3r^) r ,Di < r < D 2 
[0 ,r>D2, 

(71) 

where r is the center-center distance between 
two particles, as well as a repulsive potential. 
Use = esc(2a/r)^"‘. The values of the parame¬ 
ters in these potentials are set to Di = 2.245a, 
D 2 = 2.694a, c = 58.5/a^, and esc = 10.0. As 
the cluster evolved, they monitored the evolu¬ 
tion of the radius of gyration 


Rg{t) 


1 


ly p 

E 


1/2 




(72) 


AU/keT 

c 

esc a/Ax 

1] AtjtDa Np Ng 

2.39 

58.5/a^ 

18.0 3.29 

1 2 • lO""" 13 32^ 


Table 1. Parameter values for the icosahedron col¬ 
lapse simulations. 


where is the position of the center of mass. 
Their results revealed that the hydrodynamic 
interactions slow the time of cluster collapse and 
lead to a final state in which the particles rear¬ 
range themselves within the cluster. Delong et 
used this test case to validate the treat¬ 
ment of hydrodynamic interactions in fluctuat¬ 
ing IBM, showing that they matched Brown¬ 
ian Dynamics simulations with hydrodynamics 
given by the Rotne-Prager-Yamakawa tensor. 

We perform similar simulations of cluster col¬ 
lapse to both show that fluctuating FCM re¬ 
covers previous results, but also to quantify the 
effects of the stresslet included in fluctuating 
FCM. The simulations are performed using pe¬ 
riodic boundary conditions, applying directly 
the spatial discretization scheme described in 
Section [V A| In our simulations, we set the mag¬ 
nitude of the repulsive potential to esc = 18.0, 
a slightly higher value than the one used by Fu¬ 
rukawa and TanakeP^. The values of the other 
parameters used in our simulations are provided 
in Table HI We also note that the inclusion of 
the particle stresslets necessitates a numerical 
scheme that accounts for Brownian drift such 
as the DC or RFD. If the stresslets were not in¬ 
cluded, the equations of motion could be inte¬ 
grated using the Euler-Maruyama scheme as the 
Brownian drift term is zero^ in an unbounded 
or periodic domain. 

Fig. 0 shows the time evolution of the av¬ 
erage radius of gyration given by fluctuating 
FCM with and without the particle stresslets. 
We also show data from simulations using the 
same interparticle potentials, but with hydro- 
dynamic interactions ignored completely. For 
each case, the averages are obtained from 150 
independent simulations. In addition, when the 
stresslets are present, we integrate the equations 
of motion using both central RFD and the DC. 
As in Furukawa and TanakaP^ and Delong et 


























21 



10 '^ 10 ’ 10 ° 10 ^ 10 ° 


tjtB 

Figure 7. The radius of gyration Rg{t) (in units of Aa;) during the time interval t/toa. = 0 — 150. The 

results are obtained by averaging over 150 independent simulations. - : no hydrodynamic interactions 

(HI); - : FCM without stresslets, Euler-Maruyama scheme;.: FCM-S, central RFD; - : FCM-S, 

DC. Error bars are shown for FCM without stresslets and for FCM-S with the DC. The inset on the left 
shows the initial configuration of the particles. The dark purple sphere is initially at the center of the 
icosahedron. The insets on the right show the spheres at one realization for the simulation with no HI and 
for FCM-S with the DC, respectively. 


we see from Fig. [^that the hydrodynamic 
interactions impact both the evolution and fi¬ 
nal value of Rg by slowing down the collapse of 
the cluster. Note that in the limit < —>■ oo, all 
the simulations must go to the same final Rg 
since the equilibrium state does not depend on 
the dynamics but only on the interaction po¬ 
tential. Here, we see also that the inclusion 
of the stresslets leads to a higher value of Rg 
during the collapse. The value of Rg with the 
particle stresslets given by central RFD and the 
DC are nearly identical. We note that for each 
timestep, the central RFD scheme requires the 
stresslets to be solved for three times. One cen¬ 
tral RFD timestep therefore requires approxi¬ 
mately 30 Stokes solves. While providing sim¬ 
ilar results as the central RFD approach, the 
DC requires just one iterative solve per time 
step (about 10 Stokes solves). This is confirmed 
from our simulations which show that the av¬ 
erage simulation time with the DC is approxi¬ 
mately three times less, « 3, than 

that for central RFD. Since the 6 parameter is 
only limited by roundoff errors, forward RFD 
would provide similar results while saving one 


stresslet iteration per timestep. For forward 
RFD, the simulation time would be double that 
of the DC. Again, we note that compared to a 
completely deterministic simulation, the addi¬ 
tional cost per time-step of including Brownian 
motion for FCM with the DC is just the distri¬ 
bution of the random stresses on the grid, an 
0{Ng) computation, and one additional Stokes 
solve, which for our FFT-based solver incurs an 
0{Ng log Ng) cost. 


2. Aggregation and percolation in colloidal 
suspensions 

As a final test and demonstration of the DC, 
we perform a series of fluctuating FCM simu¬ 
lations to examine colloidal gelation of a sus¬ 
pension of Brownian particles. We compare our 
results with those given by accelerated Stoke¬ 
sian Dynamics (ASD) in Cao et al^. While we 
do not incorporate near-held lubrication hydro¬ 
dynamics as is the case with ASD, we hnd the 
huctuating FCM provides a very similar charac¬ 
terization of the gelation process, and with the 
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DC we find the computations take a fraction of 
time of the ASD simulations. 

As in Cao et we employ an interparticle 
potential that is a combination of a (36 — 18) 
Lennard-Jones-like potential and a repulsive 
long-range Yukawa potential, 


Ujr) 

ksT 


= A 


36 


+B 


2a 
r 

exp [—K (r — 2a)] 


18 


(73) 


Along with the other parameter values, the ex¬ 
act values of A, B, and k used in our simulations 
are provided in Table [llj These parameter val¬ 
ues match exactly those used by Cao et 
including the very small time-step that must 
be taken due to the stiffness of the Lennard- 
Jones-like potential. To increase the volume 
fraction, (p, we increase linearly the number of 
particles, Np while keeping the computational 
domain over which we solve the Stokes equa¬ 
tions a fixed size such that p = N^Aira^ /(3L^). 
This is different from Cao et where the 
number of particles was kept fixed as the vol¬ 
ume fraction was varied. 

As each simulation progresses, we monitor 
the number of bonds per particle, Nh/Np. We 
use the criteriorfi^that a bond between two par¬ 
ticles is formed when their center-to-center dis¬ 
tance is less than or equal to 2.21a. Fig. 8a 
shows the time evolution of N^/Np for differ¬ 
ent values of p. We observe for a fixed time, 
the number of bonds per particle increases as 
(j) increases. We do see, however that as time 
increases, all simulations approach the same 
asymptotic value of Nh/Np « 3.47. These ob¬ 
servations consistent with the ASD results ob¬ 
tained by Cao et 

The aggregation dynamics can also be quan¬ 
tified by the time evolution of the number, N^, 
of particle clusters. We determine Nc{t) by first 
compiling a list of bonded particle pairs at time 
t and then processing the list using the fc—clique 
percolation algorith m im plemented in Python 
by Reid et a/.ESl Fig. gb shows Nc as a function 
of time. We see that for each value of p, the 
number of clusters decreases with time. We also 


find that for a fixed time, the number of clusters 
decreases as (p increases. The final structures we 
observe at t/to^ = 300 for different values of p 
are shown in Table [Till For all values of p except 
p = 0.04, we find that the particles aggregate 
to form a single structure that spans the entire 
domain. For p = 0.04, which is the most dilute 
case we considered, the final state consists of 2 
clusters that are not long enough to connect the 
opposite sides of the domain. The tendency to 
form a single cluster was also noted by Cao et 
a/.E^ and Furukawa and TanakrP^. They also 
showed that for low volume fractions, hydro- 
dynamic interactions between the particles are 
needed in order to capture the formation of a 
single percolated structure. 

In addition to achieving very similar results 
as ASD, fluctuating FCM with the DC integra¬ 
tion scheme provides these results at a relatively 
low computational cost. The average compu¬ 
tational time needed for fluctuating FCM with 
the stresslets and the DC to reach t/tn^ = 100 
(10® time-steps) with Np = 1674 was 2.5 days. 
The ASD simulation^^^ required approximately 
10 days to reach t/to^ = 100 (10® timesteps) 
for half as many particles {Np = 874). Our 
approach is 8 times faster when run on 64 In- 
tel(r) Ivybridge 2.8 Ghz cores. To further test 
the scalability of fluctuating FCM with DC, 
we simulated a suspension of Np = 3766 par¬ 
ticles at a volume fraction of p = 0.08 using 
64 Intel(r) Ivybridge 2.8 Ghz cores. Snapshots 
from the simulation at different times are pro¬ 
vided in Fig. and we include a movie con¬ 
structed from our simulations in the supplemen¬ 
tary material. This simulation required 8 days 
to reach t/to^ = 145 (1.45 x 10® time steps). 
For this larger system size, we again find that 
the particles do eventually aggregate to form a 
single structure that percolates the entire do¬ 
main. 


VI. CONCLUSIONS 


In this study, we presented a mid-point time 
integration scheme we refer to as the drifter- 
corrector (DC) to efficiently integrate the over- 
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A/ksT B/ksT 

hi 


rj a (j) Np L 

60 20a 

4/a 

10'^ 

1 3.29Aa; 0.04- 0.12 558 - 1674 128Ax 


Table II. Simulation parameters for the colloidal gelation simulations 



(a) Number of bonds per particle Ni,/Np as a function(b) Number of clusters in the domain as a function 
of t for different values of (p. of t for different values of p. 

Figure 8. Time evolution of the aggregation process for (j> = 0.04 — 0.12 and n = 4/a. 


damped equations of motion for hydrodynami- 
cally interacting particles when their Brownian 
motion is computed using fluctuating hydrody¬ 
namics. Methods based on fluctuating hydro¬ 
dynamics such as the stochastic and fluctuat¬ 
ing IBMs and the fluctuating FCM have been 
shown to drastically reduce the cost of including 
Brownian motion in simulation by simultane¬ 
ously computing the deterministic and random 
particle velocities at each timestep. 

Like Fixman’s method and RFD-based 
schemes, our DC scheme allows for the effects 
of the Brownian drift term to be included with¬ 
out ever having to compute it directly. The 
DC, however, was designed especially for the 
case where imposed constraints on the rate-of- 
strain are used to provide a more accurate ap¬ 
proximation to the particle hydrodynamic in¬ 
teractions, as is the case with fluctuating FCM. 
Imposing these constraints incurs an additional 
cost for existing schemes. Using these schemes 


with fluctuating FCM, we showed that this con¬ 
straint would at the very least double the cost 
per timestep, which in practice meant an addi¬ 
tional 10 Stokes solves per timestep. We show 
that by using the DC with fluctuating FCM, 
this additional cost can be reduced to a single 
Stokes solve per time step. Though we have de¬ 
veloped the DC for and tested the scheme with 
fluctuating FCM, FCM’s similarity with other 
methods, such as the stochastic and fluctuating 
IBMs, suggests that the DC could be an effec¬ 
tive scheme to integrate particle positions for 
these methods as well. 

Using the DC with fluctuating FCM, we have 
provided an extensive validation of the scheme, 
showing that it reproduces the correct equi¬ 
librium particle distribution, as well as pro¬ 
vides the distribution dynamics in accordance 
with the Smoluchowski equation. We have also 
demonstrated the effectiveness of fluctuating 
FCM with the DC for colloidal suspension simu- 
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Table III. State of the suspension at tjtD^ = 300 for different values of 4>. The numbers in parentheses 
correspond to the number of particles in the domain Np. The label “P” indicates that the particles 
aggregated to form a percolated network, while “NP” indicates that they have not. The snapshots show 
the suspensions at the final time and Nc. is the number of clusters in the images (accounting for periodicity). 



(a) tltoa^ = 0. (b) t/toa. = 8. (c) tjtoa = 145. 


Figure 9. Snapshots of the aggregation process for 0 = 0.08 and Np = 3766 (see also supplementary 
material^. 


lations, examining the collapse of a small cluster 
of particles, as well as the gelation process of a 
suspension. In doing so, we were able to both 
quantify the effect of including higher-order cor¬ 
rections to the hydrodynamic interactions, as 
well as demonstrate the ability to accurately 
simulate colloidal suspensions at large-scale. In¬ 
deed, we are currently applying this approach 
with active particle mod els^ t o understand re¬ 
cent experimental result^^Hm on the dynamics 
of passive Brownian tracers in active particle 
suspensions. 


The methods described in this work may also 
be modihed and extended to simulate particles 
with more complicated shap eJ^. While one ap¬ 
proach is to construct rigicPSMl^ or flexible^ 
particles from assemblies of spherical particles, 
another approach is to modify the shape of 
the kernel in the projection and volume av¬ 
eraging operators to reflect the shape of the 
particle itself. This has been done with FCM 
for ellipsoidal particle^ and, as the underlying 
framework of the method remains unchanged, 
fluctuating FCM with the redefined operators 
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should yield particle velocities that satisfy the 
fluctuation-dissipation theorem. 

A current challenge is to adapt the DC to 
these particle shapes as the Brownian drift term 
will also depend on the orientation of the par¬ 
ticle. It would be of interest to compare fluctu¬ 
ating FCM results with the experimental mea¬ 
surements of Han et aim and the recent simu¬ 
lations of De Corato et aim. We are currently 
pursuing this line of research to provide a simi¬ 
lar set of efficient and accurate simulation tech¬ 
niques to understand the dynamics of suspen¬ 
sions of ellipsoidal Brownian particles. 
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Appendix A: Stress-free boundary conditions with 
fluctuating hydrodynamics 

We have seen in the text that the slip bound¬ 
ary conditions can be imposed by requiring that 
the fluctuating stress satisfy certain symmetry 
conditions. In this appendix, we show analyt¬ 
ically that the fluctuating stress with modified 
statistics yields flow correlations proportional to 
the appropriate Green’s function. 

For the case where the fluid is in the positive 
half-space, i.e. z > 0 , with the conditions that 
Uz = 0 and duxjdz = duyfdz = 0 at z = 0 , the 
Green’s functiorP^ for this Stokes flow is 

= Gy (x-y) -h 7 jfcGifc(x-Y) (Al) 


where 


yjk — ^jk 

Y = y-2(y-i)i. 


Gy(x) 


1 / ^ij XiXj \ 

Stt?? v|x| |x|3y' 


(A2) 

(A3) 

(A4) 


Thus, if there is a point force, F = {F^, Fy,Fz), 
located at the point y in the flow domain, 
the boundary conditions can be satisfied by 
introducing an image point force F*"* = 

{Fx, Fy, —Fz) below the slip surface at Y. 

This image system approach to obtain the 
Green’s function readily extends to the case 
where there are two parallel no flux, slip sur¬ 
faces at z = 0 and z = Lz, and periodic bound¬ 
ary conditions at a:: = 0 and x = and y = 0 
and y = Ly. For our discussion, we consider the 
specific case where Lz = = Ly /2 = L/ 2 , 

but it can be easily generalized to different do¬ 
main sizes. For the case where the fluid do¬ 
main is given by z G [ 0 , L/ 2 ), x G [ 0 , L), and 
y G [0, L), the Green’s function can be ex¬ 
pressed as 





kiki 

1^ 




To impose the slip condition, we would like 
then for the fluid flow resulting from the fluctu¬ 
ating stress to have the following statistics 


(u(x)) = 0 (A 6 ) 

(u(x)u'^(y)) = 2 /cBrL(x,y). (A7) 


We show here that this is achieved by having 


(P,,(x)) = 0 (A 8 ) 

(Py(x)Pfe/(y)) = 2kBTAijkiS{yi - y) 

+2fcBTr,,w(5(x-Y)(A9) 


where 


AijJil — ^ik^jl F ^il^jk (-^ 40 ) 

Lijki = liklji + liiljk- (All) 

We notice that the fluctuating stress is L- 
periodic in each of the three dimensions. We 
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can therefore use Fourier series to obtain the 
flow resulting from the fluctuating stress. We 
define the Fourier transform as 

g(x)=^g(k)e*'^" (A12) 

k 


and, consequently, 


g(k) = L-3 / (A13) 

Jn 

The flow in Fourier space is given by 

^ (A14) 

and, therefore, 

(u,(k)u„(q)) = -0^ - 0^ 

X 0-^-0 {PA^)Ppr{^))■ 

(A15) 


From the definition of the Fourier transform and 
Eq. (A91, we find that 


(Pij(k)Pfei(q)) = 2kBTL 

+2kBT L~^T 

(A16) 


where K = — 7 k. Using this expression in (A151 
and the fact that = fc^, we find 


(ui(k)u„(q)) = 


2fcBr 

r]L^k'^ 

2kBT 

rjL^k'^ 


din 


hkr, 


+ 


lin + 


kiK, 




(5q,-k 
<5q,K- 

(A17) 


In real space, the flow correlations are then 


({ti(x)M„(y)) = 




ik x iq v 
ke 


+ 


kiKn 


fc2 




(A18) 


Performing the sum over q for each term, we 
have 


(Mi(x)u„(y)) = 



({ij(x)'u„(y)) = 


E 

k^O 


2kBT 

rjL^k^ 


_ kjkn \ ,k (x-y) 
Ozn ^2 ' ^ 


-E 

k^O 


2kBT 

riL^k^ 


din + 


ki kn 

~w 


line 


,ik.(x-Y) 


(A19) 


Writing this expression slightly differently as 
2kBT 1 / kiki 


{Uz{x.)Un(y)) = 


PS 2-^ 
k^O ' 

„-ik y ^ 




we immediately see that 

(u(x)u^(y)) = 2fcBTL(x,y), 


(A20) 

(A21) 


with L given by Eq. (A5), as desired. 


Appendix B: Error expansion for the DC 


In this appendix, we present the error analy¬ 
sis showing that the DC is weakly accurate to 
first-order in time. To simplify the analysis and 
make the presentation more compact, we write 
the scheme using the mobility matrices as 


yk+1/2 ^yk ^iik 


yk+l 


(Bl) 

-pfc+1/2 


(B2) 
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where 


i)fc+l/2 _ ^fc+l/2|-Qfej 

cfc+1/2 'T^£S\k-\-l/2pk-\-\l2 

^ — '^FCM ^ ’ 

£fc+l/2 = _/C'=+l/2[ufc], 


(B3) 

(B4) 

(B5) 

(B6) 


and satisfies Eq. ( [2^ for the realization of 
the fluctuating stress at time tk- We now ex¬ 
pand the quantities on the right-hand side of 
Eq. (B2) about t = tk, and using Eq. (Bl) we 
obtain 


V^+l/2 = Vfe + 


■pk+1/2 _ -pk 

r\ *' m 


2 dyp 
At dTl 


n'.kn'.k 

8 ^ 


2 dyp^^^ 8 dy^dyp^^^'^ 


0(At) 

-0(Af3/2) 


.,vj^-,k+i/2_ ,,vF-,k -. 

■^a/3 +T dy ^7 


8 




0(At3/2). 


(B7) 

(B8) 

(B9) 


where we have used Greek letters 
vectors and mobility matrices and 
tion that their repetition implies 
The expansions for T'^'+i/^ and Ai 
identical to those for J^'^+i/2 and 
respectively. For clarity, we have 
label '"FCM — S” for the mobility 


to index the 
the conven- 

summation. 
Vr;fe+l/2 


FCM-S 


are 


FCM-S ’ 
dropped the 
matrices in¬ 


volved in the expansion. We also have from Eq. 


(47) that 


v'^ = 


AtdUl 

^Wo.' 


(BIO) 


Inserting the expansions into Eq. (B2), the ex¬ 
pansion for the increment AJ^^ = 3^+1 — 


IS 


Aj;^ = AtV^ + At 

+Ae 


+ ^( + V." 


dyp 



1 d ( ^,VF-,k-rk I f.kfik 


+ 0{At^). (Bll) 
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Taking the ensemble average of this expression 
and using the fact that 


d 

Wp 


(yluf) 


dV^ 

dyp 


ul + v, 


ai 


we have 

{Ay^j = At 


^VF-,k pk 
■^aP -^P 


M 




kq-k 

Ol^ ' 0 

pO{At^). 

(B13) 
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Since 

(Va^p) = {^a^p) ~ ■^FCM-,a'y'^FCM-~/si^S^p)’ 

(B14) 

with the correlations on the right-hand side 
giverPSI by 

= ‘^-^MYcI,cP (B15) 

(Y^p) = -^FCM-SP (B16) 

we have that 

{YUp) = (B17) 


We emphasize that the mobility matrix in Eq. 
(B17) is the stresslet-corrected mobility matrix, 
the expression for which is given in Eq. (27). As 
a result, we obtain the correct first moment up 
to first order in time 


{Ay^j = At 


^VF-,k -pk 
^aP -^P 


A/fVT ■,kq-k 

■^aP 'P 


+AtfcBr Y +o(At"). 

oyp 

(B18) 

The outer product of the increment gives 


1. Simplified scheme for periodic, no-slip, and 
no shear boundary conditions 


From the analysis above, we see that the in¬ 
clusion of is needed to introduce 'sly -U in 
order for the scheme to account for the diver¬ 
gence of the mobility matrix. If, however, we 
have that Vy -U = 0, then we may set = 0. 
We show here that this is the case when peri¬ 
odic, no-slip, or no flux boundary conditions are 
imposed at the fluid boundary. We begin by 
noticing that the contribution to Vy • U from 
particle n is 


dY 


\ 9A(x-Y")^3 
dY- ^ " 
9A(x-Y")^3 


dxi 


d-^x 


- j ^ (uiA(x-Y”))d^x 

-b /^A(x-Y")d3x. (B22) 


where repeated Latin indices imply summation. 
The second integral on the right hand side of Eq. 
( B22| ) is zero since V • u = 0. After applying 
the divergence theorem to the first integral on 
the right-hand side, we have 


+AtYYYp + vYi 


o{Ae) 

(B19) 


where 


nk \ yiS^SF-.k-i~k I \ ^VF^kyj-k , ^ {'\lkFjk\ 

ya = M^p J^p+M^p 

(B20) 

Taking the ensemble average of the outer prod¬ 
uct gives 

{Ay’iAy'D = Ae{vYp) + o(A<") 

= 2kBTAtMlY 

+0[AY, (B21) 

showing that the second moment is also recov¬ 
ered to first order in time. 


= - y fi*n,A(x - Y")d5, (B23) 

where rii is the normal pointing out of the do¬ 
main. Thus, if UiUi = 0 pointwise, which is 
the case for no-slip and no-flux boundaries, or if 
periodic boundary conditions are imposed, then 
this integral is also zero. As this will be the case 
for all particles, we then have that 

VyU = Q. (B24) 

Thus, as a result, we may set v'^ = 0 and the 
DC scheme will still account for the Brownian 
drift. 

As this result relies on using the continous 
spatial operators, it is important to check that 
it also holds when the Stokes equations and the 
projection and volume averaging operators are 
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Figure 10. (|V ■ U{z)\)yfKi along a in a periodic 
domain and between two shear free surfaces. The 
expectation is taken using 200 realizations of the 
flow. 


discretized in space. In particular, we need to 
ensure that dUp / dY/^ is nearly zero in the sim¬ 


ulations. Fig. 10 shows {\dUp/dYp\) as a func¬ 


tion of the z-coordinate for a single particle in a 
periodic domain and between two no-shear sur¬ 
faces. For the periodic case, we see that this 
quantity is 0(10“®) and does not depend on 
z. We may therefore safely set = 0 in our 
periodic simulations. For the no-shear bound¬ 
aries, however, we observe a clear dependence 


J 


of {\dUp/dYpl) on z near the boundaries. Fur¬ 
ther investigation revealed that this dependence 
is linked to the sum used to approximate the in¬ 
tegral in Eq. (481 as we found that {\dljp/dY^\) 
is well approximated by the error estimate for 
the quadrature scheme (not shown). Thus, in 
our simulations with no-shear boundaries, we 
retain as part of the time integration. We 
note, however, that for the computation of 
one only needs to consider the particles that 
overlap the no-shear boundaries. 

Appendix C: Energy balance for the 
force-coupling method 


For a suspension of rigid particles in Stokes 
flow, the rate of work done by the particles on 
the fluid is balanced by the viscous dissipation 
in the fluid and the rate of work done by the 
external boundary. Thus, 



e : ed®x — 


• cr • MS = ^ F" • V’" 

n 

n 

(Cl) 


where n is the unit normal to the external 
boundary. For FCM, we would like to ensure 
that our definitions of particle motion yield the 
same energy balance. Taking the dot product 
between u and Eq. ([^ without the fluctuating 
stress, and integrating over the fluid domain, we 
find that 


J e : ed®x — J u • cr • MS = ^ F" • J uA„(x)d®x -|- ^ ^ r" • J ^ ^ V0„(x)d®x 

+ / (uV0„(x)-f-(uV0„(x))^) d®x (C2) 

n 


where we have also integrated by parts to obtain 
the viscous dissipation and boundary term on 
the left-hand side. 


We see that in order for Eq. (C21 to com ply 
with the standard energy balance, Eq. (Cl I, we 
need to require that the particle velocities are 
given by Eq. 0 and the angular velocities are 
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provided by Eq. (18). Additionally, in order for 
the stresslets to perform no work on the fluid, 
the local rate of strain given by Eq. (191 must 
be constrained to be equal to zero. 


Appendix D: Interpolation of the 
stresslet-corrected mobility coefficients 


In order to rapidly generate a large number 
of statistics to obtain the distribution of a sin¬ 
gle particle in a slip channel, we determine the 
motion of the particles using the FCM mobility 
matrices for the single particle and the volume 
averages of the fluctuating fluid flow. This al¬ 
lows the usage of a single realization of the fluc¬ 
tuating fluid flow to determine the velocity of 
many non-interacting particles. 

The velocity of particle n is given by 

V„ = (Dl) 


where is given by Eq. ( [^ with the 

stresslet corrected coefficients, and 

the matrix relates the particle velocity 

and local rate of strain. The quantities V„ and 
En are the random velocities and local rates of 
strain obtain by applying the FCM volume av¬ 
eraging operators to the fluctuating flow field. 

Due to symmetry, is diagonal (see Eq. 

(^1), and only and = Cy fy are non¬ 
zero. The values, however, for these non-zero 
entries depend on the particle’s distance from 
the boundaries. Thus, to perform the simula¬ 
tions, we compute the non-zero matrix entries 
at 400 equi-spaced positions between z = a and 
z = Lz — a and use linear interpolation to find 
the values at an arbitrary position along the 
channel. 
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